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Using electronic structure calculations, we conduct an extensive investigation into the Hf-Ta-C system, which 
includes the compounds that have the highest melting points known to date. We identify three major chemical 
factors that contribute to the high melting temperatures. Based on these factors, we propose a class of materials 
that may possess even higher melting temperatures and explore it via efficient ab initio molecular dynamics 
calculations in order to identify the composition maximizing the melting point. This study demonstrates the 
feasibility of automated and high-throughput materials screening and discovery via ab initio calculations for 
the optimization of “higher-level” properties, such as melting points, whose determination requires extensive 


sampling of atomic configuration space. 
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High-performance refractory materials [1-5] play an im- 
portant role in applications ranging from gas turbines to 
heat shields for hypersonic vehicles. With melting points 
above 4000 K, hafnium carbide [6-12] and tantalum car- 
bide [6,13,14] are among the most refractory binary com- 
pounds known to date [15]. Their mixture Ta,HfCs melts 
at 4215 K [16], which has long been considered the highest 
melting temperature for any solid [17]. Very few measurements 
have been documented, because of the obvious experimental 
difficulties at extreme temperatures. 

Computational approaches to melting point prediction 
offer exceptional control and monitoring of thermodynamic 
variables [18,19] and can more flexibly handle a wide range of 
materials and temperatures. However, a melting temperature 
calculation from quantum-mechanical methods has long been 
considered a challenging task, due to the requirement of 
extensive sampling of atomic configuration space, particularly 
for the liquid phase [20-22]. Given the computational burden 
of DFT [23], it is extremely difficult to perform systematic 
and high-throughput melting temperature calculations directly 
from first principles. We recently developed the small-size 
coexistence method [24,25], and managed to reduce the 
computer cost drastically. An automated computer code is 
prepared and freely distributed for direct DFT melting point 
calculations [25]. In this work, we apply the code to study 
these most refractory materials. We demonstrate that it is 
feasible to perform high-throughput materials screening 
and discovery directly via ab initio calculations for the 
optimization of melting point. 

To help identify the factors leading to high melting points 
and validate our computational methodology, we first explore 
trends among known classes of refractories. Our investigation 
focuses on the rocksalt structure in the Hf-Ta-C systems, 
because it is the only stable solid-state form in the temperature 
region of melting [26,27]. Employing the small-cell coexis- 
tence method [24], we calculate the melting temperatures of 
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rocksalt HfC, (x € [0.75,1]), as shown in Fig. 1. Our calcu- 
lations successfully capture the volcano-shape melting curve, 
which is widely observed in experiment [7,9—1 1], as well as the 
location of the apex (within 3 at. % of carbon content). Starting 
from stoichiometric HfC, the melting temperature increases 
along with carbon deficiency, until it reaches a maximum at 
the congruent melting point, near 45 at. % C (HfCo,.32). This 
feature explains why HfC undergoes carbon loss when it is 
heated and melted [29]: it remains in the solid state and loses 
carbon as the temperature increases. Further decrease in carbon 
composition leads to a drop of melting temperature. 

It should be noted that our calculations identify the so-called 
midrib surface, that is, the temperature where the free energy 
of the solid and the liquid are equal at a given composition. 
The midrib surface generally lies between the solidus and 
liquidus, but agrees exactly with the melting point at extrema 
(where the solidus and liquidus meet), thus enabling us to 
accurately predict optimal melting points. As the calculations 
are performed under constant pressure conditions, they can 
determine whether the solid melts or sublimates and only 
melting is observed. However, the calculations do not include 
the possible effects of an oxygen-rich environment (carbon loss 
via oxidation) and are thus representative of heating under an 
inert atmosphere (e.g., nitrogen). 

We note that the vertical shift of the calculated melt- 
ing curve, relative to experiment, does not appear to be 
composition-dependent and thus does not significantly affect 
trends. This shift amounts to about 5% of the melting temper- 
ature itself, which is typical for DFT calculations. In addition 
to the Perdew-Burke-Ernzerhof (PBE) functional [30], which 
we employed for melting temperature calculations, we have 
cross-checked a subset of data points with a generally more 
accurate, but considerably more expensive, Heyd-Scuseria- 
Ernzerhof (HSE) hybrid functional [31] and found an average 
shift upward by +460 K (see Table I in Ref. [25]). The fact 
that PBE and HSE-based results bracket the experimental 
data further verifies that the error is mostly caused by the 
drawback of DFT exchange-correlation functionals, which is 
understandable given the nontrivial electronic structure of HfC 
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FIG. 1. (Color) The Hf-C phase diagram. Prior experimental 
measurements of the melting points are compared with the present 
computational results (labelled “PBE” as they rely on the PBE 
functional). The temperatures where free energies of the liquid 
and of the solid intersect are marked by “+”. Also shown are the 
solidus and liquidus obtained via a CALPHAD model [28] fitted to 
our calculated thermodynamic data. Our calculations successfully 
capture the location of the apex within 3 at. % of carbon content. The 
vertical shift in calculated temperature, relative to experimental data, 
is essentially composition-independent and mostly reflects DFT error 
(see discussion in text). 


and the extremely high temperature. Despite the considerable 
difference between PBE and HSE results, we note that all 
melting point calculations are based on the PBE functional in 
this article. Therefore the trend of melting temperature change 
is still valid and consistent, and the comparison of melting 
temperature remains effective, as all calculations are carried 
out with the same DFT functional. 

The high melting temperature of hafnium carbide is pri- 
marily (through the well-known relation T„, = AH/AS) due 
to its exceptionally large fusion enthalpy of 0.81 eV atom™!, a 
value usually unparalleled among refractories. (For reference, 
Al,O3 (m.p. 2345 K): 0.22; W (m.p. 3695 K): 0.37; Hf 
(m.p. 2506 K): 0.26 eV atom™!). Indeed, a large heat of 
fusion is the first and most prominent factor we find that 
contributes to a high melting point. The chemical origin of 
the remarkably large heat of fusion can be studied via a 
wave-function analysis. While most researchers agree that the 
Hf-C interatomic bonding is a mixture of metallic, covalent 
and ionic interactions, its precise nature has not been well 
understood [32-34]. The system’s wave functions, illustrated 
in Fig. 2, reveal numerous types of chemical interactions, 
including Hf—Hf 5d o bond, C-C 2p o bond and Hf 5d-C 2p 
x bond. This diversity enables each atom to bind with all its 
first- and second-nearest neighbors, thus forming an unusually 
large number of bonds and promoting the formation of a 
deep valence band (see Fig. 2). These bonds also carry both 
covalent and ionic characters. On one hand, the decomposed 
density of states (Fig. 2) shows contribution from both carbon 
and hafnium, hence demonstrating a typical covalent bond 
pattern. On the other hand, charge density analysis (Fig. 3) 
clearly shows partial charge transfer from hafnium to carbon, 
a strong evidence of ionic bonding, which is confirmed by a 
Bader charge analysis [35] indicating a 0.62 e charge transfer. 

The second contributor we recognize is the presence of 
point defects, which affect melting temperature via entropy. 


RAPID COMMUNICATIONS 


PHYSICAL REVIEW B 92, 020104(R) (2015) 


5 
© 
ð 


Density of states 


-12 -10 -8 -6 -4 -2 0 2 4 
Band energy (eV) 


FIG. 2. (Color) Electronic density of states in HfC showing clear 
participation of both Hf and C in forming covalent bonds. Total 
density of states is shown in black. Projections on C and Hf are 
colored in red and blue, respectively. The Fermi level is at 0. Vertical 
lines are energy levels of atomic orbitals. Insets are wave functions 
illustrating the diversity of bond types in HfC with clear covalent 
character. From left to right, these figures represent Hf 5d o bond, C 
2p o bond, and Hf 5d—C 2p x bond. Surfaces of constant value of 
the real parts of the wave functions are represented. Hf and C atoms 
are colored in green and blue respectively. 


More generally, we find that, at high temperatures, entropic 
effects favor and stabilize a considerable amount of lattice 
defects. When solid HfC becomes off-stoichiometric HfC,_,, 
the presence of carbon vacancies increases the configurational 
entropy (e.g., for an ideal lattice solution, S = —k a xi Inx;), 
and this benefit is further magnified by the high temperature 
(G = H — T S). If this entropic effect more than offsets the 
defect formation energy penalty, these vacancies stabilize the 
solid phase. Since, by definition, vacancies can only exist in 
the solid phase, this effect is absent in the liquid and the net 
effect would be an increase in melting temperature. While this 
argument appears contradictory to the Lindemann melting 
criterion [36], we note that it is actually complementary. By 
empirically correlating melting with the amplitude of thermal 
vibration, the Lindemann’s rule focuses on the solid structure, 
and it lacks a thermodynamic awareness of both the solid and 


FIG. 3. (Color) The electron transfer in HfC. The figure reports 
p/po — 1, where p is electronic charge density from DFT wave 
function and po is initial overlapping atomic charge density. The 
sharp contrast clearly shows a charge transfer from Hf to C indicative 
of ionic interactions, although covalent character is still visible in the 
anisotropy of the charge density. 
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the liquid. Moreover, an overall evaluation needs to account 
for both energetic and entropic factors, in addition to the 
geometric structure. Though vacancies may facilitate melting 
by allowing larger displacements according to the Lindemann 
criterion, our calculation suggests that vacancies also increase 
the heat of fusion of HfC at low vacancy concentrations (see 
Fig. 2 in Ref. [25]). Therefore both energetic and entropic 
effects are in favor of a more stable solid phase in the presence 
of vacancies through a mechanism that acts independently of 
Lindemann’s suggestion. 

Indeed, this fact explains why melting point climbs 
when HfC becomes off-stoichiometric and carbon-deficient, 
a phenomenon widely observed in experiments (Fig. 1). 
Furthermore, this entropic effect becomes so large at high 
temperatures that it not only stabilizes defects, but facilitates 
their formation as well. For instance, we observe the formation 
of Cz (two carbon atoms near one anion lattice site) and va- 
cancy in MD simulations, especially for compositions close to 
stoichiometric HfC. These unstable C2 complexes tend to leave 
the solid, which results in carbon loss of stoichiometric HfC. 

The third chemical factor we identify is well exemplified 
in the Hf-Ta-C system. While binary carbides, such as HfC 
and TaC, are constrained by the given electronic properties 
of these metals, mixing two carbides provides an avenue to 
tune chemical properties. Hf and Ta share a similar electronic 
structure but a slightly different number of valence electrons, 
which allows tuning of the location of the Fermi level so that 
it lies precisely between the bonding and anti-bonding bands 
without distorting the density of states. 

To investigate this effect, we calculate melting points of 
Hf,.Taj_,Co.g75 at various compositions (x = 0,0.125,0.25, 
0.375,0.5,0.75,0.81, 1). Our choice of carbon content is based 
on a consensus of experimental data that maximize the 
melting points of the binary carbides. Because of the rapid 
evaporation of carbon at high temperatures, it is difficult to 
assign the measured melting point to the correct composition 
in experiments [12,29]. The maximal melting temperature in 
the Hf-C system falls in the HfCo.85—0.95 region [7,9, 10, 12,29], 
while Ta-C has a strong tendency to lose carbon and the 
melting point maximum is near TaCo-0.9 [12,14,29]. The 
Hf-Ta-C system contains HfTa,Cs, which has long been 
considered the most refractory substance known to date [17]. 
Our computational results, in general, agree with Agte’s 
experimental measurements [16], as shown in Fig. 4. Our 
calculation indeed captures a cusp near HfTa,Cs in the 
composition-dependence of the melting point. Our calculation 
provides more details at the Hf-rich region, which was not 
explored in Agte’s experiment. Existing experimental data are 
controversial on the issue of whether HfC or TaC has a higher 
melting point. Agte found HfC and TaC have nearly the same 
melting temperature [16]. While Rudy’s study showed that 
the melting point of TaC is about 50 K higher than that of 
HfC [14], Emeleus, on the contrary, claimed that HfC melts 
at least 100 K above TaC [15]. Our calculations suggest that 
the melting temperature of HfCo.375 is ~70 K higher than 
that of TaCo.375. Given the difficulty in measuring melting 
points at such high temperatures and the melting point’s 
sensitivity to the exact carbon content, which is also hard 
to control at high temperatures, these discrepancies are not 
surprising. 
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FIG. 4. (Color) Melting temperature of Ta, Hf,_,Co.g75 as a func- 
tion of x. Agte’s measurements [16] focused on the Ta-rich side, while 
the rest part of the phase diagram was extrapolated as dash lines. 
Our melting point calculations provide more details at the Hf-rich 
region. The calculated melting curve captures the cusp near TayHfCs 
(as the melting points of Tag. 75Hfo 25Co.g75 and Tagg: Hfo,19Co.g75 are 
evidently higher than those of Tao sHfo.sCo.s875s and TaCo.g75), which 
illustrates the effect of tuning the Fermi level via alloying, so that 
it lies precisely between the bonding and antibonding bands. The 
inset shows the effect on the Fermi level of the solid phase by 
tuning composition in Ta, Hf,_,C. Vertical lines are Fermi levels. Our 
calculation suggests that the melting temperature of HfTayCs, instead 
of being the highest, barely falls within the melting temperatures of 
the pure components, which corroborates Rudy’s report [37]. 


Guided by the three contributing factors discussed above, 
we explore a new class of refractory materials, which may 
have higher melting temperatures than Hf-Ta-C. (1) We focus 
on isostructural alternatives to Hf-Ta-C, because the strong 
binding with both first and second-nearest neighbors is a 
very favorable feature we wish to maintain. For the same 
reason, we look for an alternate composition with similar 
cation/anion atomic radius ratios and similar electronegativity 
differences. (2) To preserve the ability to tune the Fermi 
level, we consider more than one element on both the cation 
and the anion sublattices and start with a composition space 
including Ta, Hf, B, C, and N. Adding more transition metals 
did not appear beneficial since the valence electron density 
of Ta and Hf already brackets the optimal Fermi level for 
the rocksalt crystal structure considered here. For the anions, 
atom size and electronegativity considerations (to preserve 
the rocksalt crystal structure) lead us to limit ourselves to 
2p elements. (3) According to known melting temperatures 
of binary compounds (HfB: 2280 (decompose); HfN: 3660; 
TaB: 3360; TaN: 3370 K), we identified the Hf-C-N ternary 
subsystem as a promising candidate. HfN, a solid in rocksalt 
structure like HfC, has the highest melting temperature of all 
known metal nitrides. The high stability of both HfN and HfC 
suggests a thermally stable ternary system. Indeed, we find that 
the Hf-C-N system generally has a larger heat of fusion than 
HfC does. The Hf-C-N system was found thermodynamically 
stable [38], and the HfC-HfN mixture features complete solid 
solubility [39,40]. (4) To exploit possible entropic effects, we 
allow for vacancies on the anion sublattice. 

Our calculations indicate that the Hf-C-N system includes 
materials that have higher melting points than any other 
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FIG. 5. (Color) Melting temperatures of Ta-Hf-C-N alloys. Filled 
circles mark the calculated melting temperatures in the Hf-C and 
Hf-C-N systems while open circles show data from the Ta-Hf-C 
system for comparison. The melting temperature surface Tp (xy, xc) 
(shown as contour lines) was obtained via a regression analysis of 
the calculated melting temperatures based on a quadratic function of 
composition. See Table II in Ref. [25] for melting point data. A 2D 
version of the melting point surface is available in Fig. 1 in Ref. [25]. 


substances known to date. As shown in Fig. 5, we find a 
large number of Hf-C-N mixtures, whose melting temperatures 
are significantly higher than the Hf-C and Hf-Ta-C systems. 
These new refractory materials increase the melting tempera- 
ture record by up to 200 K. A regression analysis of our melting 
point data indicates that the highest melting point is located in 
the vicinity of xy = 0.20 and xc = 0.27 (where x subscript 
E denotes the atomic fraction of element £). We find that Ta 
does not help increase the melting temperature further. 

In investigating this broader class of systems, we have 
observed another, independent, melting point-enhancing 
mechanism. We find that the addition of nitrogen remarkably 
changes the liquid structure and renders the phase less stable, 
which hinders melting. We explain this effect as follows. A liq- 
uid is more stable at a high temperature because it can access a 
much larger phase space, which contributes to a larger entropy 
that offsets its higher energy. In particular, the liquid allows 
for a richer variety of pairwise correlations. For instance, 
while there are only Hf-X (X = C, N) nearest neighbors in 
solid-state Hf-C-N, additional pairs such as Hf—Hf and X 1-X2 
(X1,X2 = C, N) are allowed in the liquid. This is an important 
entropic benefit in favor of the liquid phase, provided these 
new pairs do not entail too much energy penalty. We find that 
the main impact of the additional nitrogen is via the unstable 
C-N and N-N pairs, which is made clear in the following two 
analyses. First, we calculate the defect formation energy of 
X-X (X =C,N) in the matrix of solid-state HfC as 


AE = E(a X-X pair on one anion lattice site) 
+ E(vacancy) — 2E(X). 


This quantity measures the energy cost to move a X atom from 
an anion sublattice site (leaving a vacancy at the site) to another 
anion (creating a X—X pair at the site). We find N-N has a much 
higher defect formation energy than C-C (5.8 versus 3.6 eV), 
which suggests a larger energy penalty when breaking Hf—N 
bonds to form N-N, a necessary step to melt the solid. As 
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FIG. 6. (Color) Pair correlation function (normalized as r —> oo) 
in liquid-state Hf32C24N7. We find a remarkable difference between 
carbon and nitrogen: there is nearly no N—N (yellow) “near neighbors” 
(first pair correlation peak), compared to a considerable amount of 
C-C (red). The inset shows compositions of near neighbors for C and 
N atoms in liquid-state Hf3.C.4N7 (the compositions of Hf—C and 
Hf-N account for the rest and are omitted). N atom has significantly 
less tendency to couple with C (2.17% vs 7.15%) and N (0.02% vs 
0.59%). 


this process becomes less favorable with nitrogen added, the 
Hf-C-N system is harder to melt. Indeed, the heat of fusion 
is larger in the Hf-C-N system (Fig. 2 in Ref. [25]). Second, 
the pair-correlation function in Hf-C-N liquid (Fig. 6) shows 
dramatically lower occurrence of C—-N and N-N, compared to 
a considerable amount of C-C pairs. This is also due to the 
higher formation energy of these two pairs. Indeed, a nitrogen 
atom has significantly less tendency than carbon (Fig. 6) to 
couple with the anions (C and N). The addition of nitrogen 
atoms largely reduces the number of anion-anion pairs in the 
liquid, forcing them to bind with Hf. This constraint limits 
the accessible phase space of the liquid and thus reduces its 
entropy. Figure 3 in Ref. [25] formally quantifies this effect 
in the Hf-C-N system. The relative instability of the N-N 
bond may appear strange at first, as the nitrogen molecule 
(N2) is usually considered very stable. However, when a N-N 
complex is bound to other atoms, its stability can significantly 
decrease, as commonly observed in other compounds. For 
example, azide, a common compound with an anion N,, is 
usually explosive; hydrazine (N2H4) is used as high-energy 
rocket fuels; dinitrogen tetroxide (N20) easily undergoes 
decomposition to its monomer (NO?) under room temperature. 

In summary, we have performed ab initio electronic struc- 
ture calculations to study the Hf-C and Hf-Ta-C systems, which 
hold the highest melting temperatures known to date. The 
ab initio methods enable us to freely explore new compositions 
in a complex multicomponent system without the demanding 
prerequisite to develop new empirical potentials. The first- 
principles electronic structure theory also naturally captures 
the electronic effects driven by the position of Fermi level, 
while empirical potentials usually fail to do so. We have 
identified and investigated three factors responsible for the 
exceptionally high melting points in a class of transition metal 
carbides: (i) the presence of a large number of strong bonds 
between both nearest and second-nearest neighbors that exhibit 
a mixture of strong covalent and strong ionic characters; 
(ii) the entropy contribution of point defects that can exist 
in the solid but in the liquid (such as vacancies); and (iii) the 
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ability to tune, via alloying, the position of the Fermi level so 
that it lies just between the bonding and antibonding bands. 
These observations suggest the exploration of the Ta-Hf-C-N 
system in order to further increase the melting point. Our 
calculations suggest that a Hf-C-N alloy with 20 at. % of N 
and 27 at. % of C increases the current melting point record by 
up to 200 K and identify a melting point increase mechanism 
mediated by changes in pair correlation functions. 
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